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We investigate ground state properties of the half-filled staggered-fiux Hubbard model on a square 
lattice. Energy gaps to charge and spin excitations and magnetic as well as dimer orders are 
calculated as a function of interaction strength U /t by means of constrained-path quantum Monte 
Carlo method. It is found that the system is a semi-metal at U/t< 5.6 and a Mott insulator with 
long-range antiferromagnetic order at U/t>^ 6.6. In the range 5.6 < U/t < 6.6, the ground state 
is an correlated insulator where both magnetic and dimer orders are absent. Furthermore, spin 
excitation in the intermediate phase appears to be gapless, and the measured spin-spin correlation 
function exhibits power-law decaying behavior. The data suggest that the non-magnetic ground 
state is a possible candidate for the putative algebraic spin liquid. 

PACS numbers: 71.10.Fd,02.70.Ss 



At sufficiently low temperatures, condensed matter 
systems have a tendency to undergo phase transitions 
and develop long range order which reflects broken 
symmetry [1]. In a two-dimensional antiferromagnct, 
however, Anderson recognized that the system could have 
a ground state that avoids all spontaneous symmetry- 
breaking and does not have magnetic order even at zero 
temperature [2]. Anderson's discovery, in conjunction 
with many subsequent theoretical investigations, uncov- 
ered a new class of matter, named spin liquids, that go 
beyond Landau's paradigm. Most notably, in contrast 
to conventional symmetry-breaking, spin liquids possess 
topological orders that cannot be characterized by local 
order parameters and carry fractionalized excitations [3]. 

Model Hamiltonians have played an important role in 
realizing such exotic spin liquid states [4, 5]. Evidence 
of spin liquid phases has been found in the spin 1/2 
Heisenberg model on triangular lattices [6], square lat- 
tices with frustration [7-9], and kagome lattices[10]. In 
these geometrically frustrated systems[ll], antiferromag- 
netic (AF) orders are suppressed by strong quantum 
fluctuations. In addition to spin systems, there is also 
progress using the Hubbard model which contains spin 
and charge degrees of freedom. Spin liquid ground states 
have been identified in the model on anisotropic triangu- 
lar lattices[12] and on bipartite honeycomb lattices[13]. 

In this paper, we examine ground state properties of 
the half-filled staggered-flux Hubbard model (sfHM) on a 
square lattice. As will be seen later, low energy physics in 
the sfHM is described by Dirac fermions, similar to that 
found in the Hubbard model on honeycomb lattices[13]. 
The model is defined by the Hamiltonian 
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where = t e is the nearest-neighbor hopping and we 
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FIG. 1. (a) Arrangement of hopping amplitudes on a square 
lattice. The phase $ = tt is distributed equally so that 
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(b) Band structure of the 



tight-binding Hamiltonian. (c) Top panel: power law expo- 
nent 77 extracted from fits to staggered spin-spin correlation 
functions. Bottom panel: thermodynamic limit charge gap 
Ac and magnetic moment (see text for definition). The 
shaded region indicates the non-magnetic insulating phase. 



set t = 1 throughout this work. The operator C;^ {cia) 
creates (annihilates) an electron with spin a J, at site 
i on a lattice of size N = L x L. f7>0is the onsite 
Coulomb repulsion. We work in the canonical ensemble. 



An electron gains a phase $ = X^n ^ij when it hops 
around a plaquette of the square lattice. $ = cor- 
responds to the original Hubbard model. We focus on 
the case $ = tt in the present study. There is a gauge 
freedom in choosing ^y. Here we distribute the phase 
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$ equally over all bonds around a plaquette and ar- 
range the hoppings according to Fig. 1(a). This leads 
to a lattice with plaquettes threaded alternatively by 
flux $ and — $. At U = 0, the energy spectrum is 
Ek = ±2-y/cos2 kx + cos^ ky. The two energy bands meet 
at the Fermi surface ek = located at nodal points 
ko = (±7r/2, ±7r/2), as shown in Fig. 1(b). Close to 
the four nodal points the energy depends linearly on k, 
which is similar to the massless Dirac spectrum found on 
the honeycomb lattice. 

Our key result is that an intermediate non-magnetic 
insulating ground state is identified between the semi- 
metal phase at weak interaction strengths and the AF 
Mott insulator at strong couplings where the hopping 
terms become irrelevant. The calculated dimer correla- 
tion function shows that columnar valence bond order is 
also absent in the intermediate phase. These results seem 
to indicate that the non-magnetic insulating phase is a 
candidate for the putative algebraic spin liquid ground 
state. Therefore, our work suggests recent progress in 
optical lattice experiments [14] might provide a promis- 
ing way of simulating the model, and observing this novel 
state of matter. 

The sfHM is solved numerically by means of the 
constrained-path quantum Monte Carlo method[15]. De- 
tails of the method are described in the supplemental 
materials. We begin with the results for the charge exci- 
tation gap. In the canonical ensemble, the charge gap at 
half-filling can be defined as[16] 
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which is the energy cost of removing a pair of electrons 
from the half-filled ground state while keeping the system 
in the Sl^^ = sector. We measure Ac(i) as a function 
of interaction strength U on lattices with linear dimen- 
sion up to L = 14. As shown in the inset of Fig. 2(a), the 
charge gap increases with U on finite lattices. To pin- 
point the critical interaction strength where the system 
turns into an insulator, we extrapolate Ac(L) at fixed U 
to the thermodynamic limit L — )■ oo using second order 
polynomials in . The results, shown in the inset of 
Fig. 2(a), indicate that the system is gapped for U > 5.6. 

In addition to the charge excitation gap, AF long-range 
order is another essential feature characterizing a Mott 
insulator. To investigate whether there is any AF order 
in the ground state, we calculate the spin structure factor 
at the Neel wave vector qAF = 
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where 5*^ is the spin operator along the ^-direction {6 = 
x,y,z), and (S^Sq) is the equal-time spin-spin correla- 
tion function. Defining m?{L) = S {q^AF , L) / L"^ , a mag- 
netically ordered phase is singled by a finite m?(L) in 
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FIG. 2. (a) Extrapolation of the charge gap Ac(i). Solid 
lines represent second-order polynomial fits to the QMC data. 
Inset shows the charge gap Ac (E) as a function of U obtained 
for L = 4, 6, 8, 10, 12, 14, and extrapolated (empty circle) 
values. Lines are guides to the eyes, (b) Finite size extrap- 
olation of the spin structure factor rr?(E). Solid lines are 
second-order polynomial fits to the QMC data. Inset: m^{L) 
versus U on finite lattices and extrapolated (empty circle) 
values. Lines are guide to the eyes. 



the thermodynamic limit. The inset of Fig. 2(b) shows 
the results of m?{L) as a function of U on finite lattices. 
We use second-order polynomials in L^^ to fit the QMC 
data and extract the value of m?[L) in the L ^ oo limit. 
It can be seen from the inset of Fig. 2(b) that AF or- 
der kicks in at t/ > 6.6, below which the system is in a 
paramagnetic phase. 

Our analysis of charge gap and magnetic order above 
suggests that the ground state of the sfHM is a semi- 
metal aX U < 5.6, and becomes a Mott insulator with 
long-range AF order at U > 6.6. Therefore, unlike 
the original half-filled Hubbard model which enters the 
Mott phase at arbitrarily small U, the sfHM has a finite 
Mott transition point. This finding is consistent with a 
previous finite-temperature determinant quantum Monte 
Carlo work[17] which reports that the critical point of 
Mott transition in the staggercd-flux model lies in the 
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FIG. 3. (a) Finite size extrapolation of the spin gap As(L). 
Solid lines are second-order polynomial fits to the Monte Carlo 
data. The inset illustrates As{L) and its extrapolated value. 
Lines are guide to the eyes, (b) Long-range behavior of the 
staggered spin-spin correlation function C(r) = ( — l)'"(5'r 5*0 + 
S^S^ + SrSS) obtained on L = 16, 20, and 24 (t/ = only). 
Straight lines are representative power-law fits to the data: 
Green (dot-dashed) line: L — 24 at U = 0, black (dashed) 
line: L = 20 ai U — 4, and purple (solid) line: L = 20 at 
U = 6. 



range 4 < J7 < 8. Moreover, our results indicate that in 
the region 5.6 < C/ < 6.6 there is an intermediate phase 
that is neither a semi-metal nor a Mott insulator. 

The absence of AF order in the intermediate phase in- 
dicates that the ground state is dominated by short-range 
spin correlations. At large distances, the spin correla- 
tion function could either decay exponentially or follow a 
power-law. To study the nature of the non-magnetic in- 
sulating phase, we first calculate the spin excitation gap. 
Following Ref. [18], we write the spin gap at half- filling 
as 
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which measures the energy cost of flipping an electron 
from spin-down to spin-up. Based on confinement argu- 
ments, a gapped spin excitation implies a finite corre- 
lation length, leading to an exponentially decaying spin- 



spin correlation. On the other hand, the correlation func- 
tion would be described by a power-law if the spin excita- 
tion is gapless. We compute As{L) as a function of U on 
finite lattices. The spin gap results are shown in the inset 
of Fig. 3(a) for 5 < [/ < 7. The data at a given U is then 
extrapolated to L — > oo using a second-order polynomial 
in to extract the spin gap in the thermodynamic 
limit. Typical behavior of the fits is plotted in Fig. 3(a). 
As expected, the extrapolated spin gap remains zero in 
the gapless semi-metal phase (17 < 5.6) and in the Mott 
phase {U > 6.6) due to the presence of gapless spin wave 
excitations. More importantly, As(L) also shows gapless 
behavior in the region 5.6 ^ U ^ 6.6, implying that the 
spin-spin correlation should follow a power-law at large 
distances. 

To support this observation, we plot in Fig. 3(b) the 
staggered spin-spin correlation function along the a;-axis. 
It appears that C(r) indeed decays algebraically at large 
separations. Moreover, the correlation function decays 
more slowly with increasing U, and starts showing sat- 
uration in the Mott phase U > 6.6. In order to quan- 
tify the long-range behavior of C(r), we fit the staggered 
spin correlation function to a power law a|r|'' for |r| > 2, 
where a and 77 are two fitting parameters. At U = 0, it is 
known that C(r) decays as |r|~''[19]. This is also demon- 
strated in Fig. 3(b): the fitted exponent of C(r) for free 
fermions on a half-filled 24 x 24 is 1] = -3.95 ± 0.13. The 
exponent as a function of U extracted from several half- 
filled lattices is plotted in the top panel of Fig. 1(c). It 
can be seen from the figure that 77 immediately increases 
with U from its non-interacting value due to the effect 
of interaction. Although the data is quite scattered, the 
figure suggests that the exponent 77 increases slowly with 
U in the region 5.6 ^ U < 6.6. 

Next we consider other order parameters proposed in 
Ref. [19]. The simplest scenario is the columnar valence 
bond solid (VBS) which breaks translational symmetry. 
The VBS order can be probed by measuring the dimer 
structure factor 



i?,,(q,L) = l^e^'i-"-C|,(r), 
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where Cgg(r) is the z-componcnt equal-time dimer- 
dimer correlation function for singlet bonds along the 
(5-direction {S = x, y) 



C|5(r) — {S^^^S^S~Sq) — {S-Sq)'^ 
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In the columnar VBS state, dimers line up coherently. 
Therefore Z?55(q, L) would pick up a characteristic mo- 
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X or 



y depending on the orientation of the bonds. Indeed 
Dss{c^,L) peaks at k^a in our finite size simulations, as 
shown in the supplemental materials. To extract the 
VBS order in the thermodynamic limit, we calculate 
dgg{L) = Dss{'k.ss)/N and extrapolate to the L 00 
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FIG. 4. Thermodynamic limit extrapolation of the columnar 
VBS order In both figures, solid lines are three representative 
third-order polynomial fits to the Monte Carlo data. Insets 
shows the normalized dimer structure for single bonds in x 
or y directions on finite size lattices. Lines are guide to the 
eyes. The extrapolated value in the thermodynamic limit is 
indicated by orange (empty) circles in the insets. 

limit. As shown in Fig. 4(a) and (b), both quantities 
vanish in the thermodynamic limit, implying the absence 
of columnar VBS order in the intermediate phase. 

A commonly adopted definition of a spin liquid is that 
it is a non-magnetic Mott insulator in which neither spin 
nor lattice symmetry is broken. Based on this defini- 
tion, our numerical data presented in this work seems 
to suggest an algebraic spin liquid ground state in the 
half-filled sfHM. However, the most unambiguous evi- 
dence of a spin liquid is its fractionalized cxcitation[20]. 
Due to the nature of our method, we arc not able to 
directly measure quantum properties of excited states. 
In terms of the method, we note that although the half- 
filled staggered-flux model does not have the fermion sign 
problem, we deliberately keep the constrained-path ap- 
proximation and calculate the ground state properties at 
half-filling. Our benchmark data shows that the error ap- 
pears to be small when compared with exact answers, as 
shown by the benchmark data in supplemental materials. 



However, it is possible that the systematic error grows 
with L. A recent exact QMC method, linearized auxil- 
iary fields Monte Carlo technique, reports that the half- 
filled ground state energy at C/ = 4 is is -0.85996(5) [21] 
per site in the thermodynamic limit. Our method, af- 
ter boundary condition averaging, gives —0.8559(4) [22], 
corresponding to a 0.47% error. By including this sys- 
tematic error, we estimate that the lower critical point 
where charge gap opens would be pushed toll ^ 5.4±0.1. 
In observable results such as correlation functions, exten- 
sive tests[15, 23] show that our systematic error is small 
even at half-filling [23] and does not affect the physics of 
the numerical solutions. 

To summarize, we have studied ground state proper- 
ties in the half-filled staggered-flux Hubbard model on 
a square lattice. Charge and spin excitation gaps as 
well as spin and dimer orders are extracted by means 
of the constrained-path quantum Monte Carlo method. 
The system is found to be a semi-metal at [/ < 5.6 
and an AF Mott insulator at C/ > 6.6. In the region 
5.6 ^ U < 6.6, our data suggests that both AF and VBS 
orders are absent in the ground state. Spin excitation in 
this region is gapless, a result that is consistent with the 
calculated staggered spin-spin correlation function which 
shows power-law decaying behavior at large distances. 
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METHOD AND BENCHMARKS 

The square lattice staggered-fiux Hubbard model 
(sfHM) is solved numerically by using the ground 
state constrained-path quantum Monte Carlo (CPQMC) 
method[l]. The CPQMC method projects the many- 
body ground state wave function |\l/o) from a trial state 
assuming (^'tI^'o) 0, by successively apply- 
ing an imaginary-time propagator e 
At being the imaginary-time step 
is decomposed according to the second order Trotter- 
Suzuki formula [2]. The two-body part of the result- 
ing operator is then transformed into one-body pro- 
jectors using a spin-decomposed Hubbard- Stratonovish 
(HS) transformation [3]. Apart from systematic errors 
due to Trotter break-up, we arrive at a formally exact 
expression e~^^^ = X]{x} ^({^})^({^})' where {x} is 
a collection of N Ising-like HS variables, P({x}) is their 
probability distribution, and i?({x}) is a one-body pro- 
jector. The projection is then realized by importance- 
sampled open-ended random walks with non-orthogonal 
Slater determinants (SDs), where the projectors i?({x}) 
propagate one SD into another. 

Away from half- filling (one electron per lattice site), 
the fermion sign problem is controlled by the constrained- 
path approximation[l]. The ground state wave function 
obtained by the CPQMC method is written as = 
'}Zij,w{4')\(t)) , where |</)) are SDs sampled by the QMC, 
w{(j)) are weight factors dictated by the distribution of 
\(j)). Since the Schrodinger equation is linear, and 
— I^q) are two degenerate solutions which can be both 
sampled in a random walk. The appearance of the two 
sets with opposite signs in the Monte Carlo samples is 
the origin of the exponential sign decay. To control the 
problem, we restrict the walkers such that at each step 
of projection the condition {^t\4') > is fulfilled. 

After the random walk has equilibrated, expectation 
values can be computed from I^'q). For example, the 
ground state energy is evaluated using a mixed estimator 



{^t\H\^1) 



(1) 



where H is the Hamiltonian. For observables that does 
not commute with H, we use a scheme called back- 
propagation (BP)[1] which is similar to the forward- 
walking technique in the Green's Function Monte Carlo 
method [4]. Throughout this work, we use free electron 
wave functions as our \^t)- A very small amount of 



anisotropy ((5$ = 0.002$) is added to the total flux per 
plaquette $ = tt in order to lift the two-fold degeneracy 
in the finite size single-particle spectrum. We also aver- 
age our results over boundary conditions so that finite 
size effects can be minimized. 

To access the accuracy of our method, we compare the 
energy of a 4 x 4 Hubbard cluster at C/ = 4 and $ = 0. 
Let Eg{Nf, Ni) denote the energy of the system with 
(Ni) spin-up (-down) electrons. The exact energy per 
site is Eg{7, 7)/N = -0.9831 [5], where TV = L x L is the 
number of lattice sites of an L x L square lattice. After 
correcting the Trotter error, the CPQMC method gives 
—0.9821(6), corresponding to an error of 0.1%. At half- 
filling, the exact energy is Eg{8, 8)/N = -0.8513[5]. The 
CPQMC energy obtained with contrained-path approxi- 
mation is —0.8491(2), which is within 0.26% of the exact 
energy. Comparisons of other observables, such as spin- 
spin correlation function, computed using the BP scheme 
can be found in Ref. 1, and 6. 



DIMER CORRELATION AND STRUCTURE 
FACTOR 

The dimer-dimer correlation function examined in this 
work is defined as 



(2) 



where is the z-component spin 1/2 operator at site 
r. d = X or y denotes the orientation of singlet bonds. 
Typical behavior of the dimer-dimer correlation function 
is shown in Fig. 1. for a 16 x 16 lattice at half-filling with 
interaction strength U — 6.8. 

We are also interested in the dimer structure factor. 
This quantity is defined as 



(3) 



Fig. 2 illustrates the dimer structure factor obtained for 
a L = 16 lattice &iU = 6.8. In this example, Dxx{^.^ 16) 
peaks at the dimer characteristic wave vector q = (tt, 0), 
indicating a weak dimer order on a finite lattice. To 
determine whether the ground state is dimerized, we ex- 
trapolate Dss{<\^ L) at q = (tt, 0) or (0, tt) to the L oo 
limit, and extract its bulk value. The results are shown 
in the main text. 
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FIG. 1. Dimer-dimer correlation function C|a;(r) obtained 

on a 16 X 16 lattice at f/ = 6.8. The correlation function FIG. 2. Dimer structure factor D^^{q,L) obtained on a L = 

is plotted along the x-axis. The inset shows the behavior of 16 square lattice at half-filling and U = 6.8. 
CLir) for |r| >2. 
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